Spatial variation of the laser fields and electron dynamics at a gas-solid interface 
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In the long wavelength domain, typically for wavelengths A ^ 100 A, the laser fields are usually 
taken as independent of the spatial coordinate. However, at the gas-solid interface the electron 
density of the material and the incident laser fields vary sharply on a scale of few angstroms. Instead 
of solving Maxwell equations, we present here a theoretical model, called Electromagnetic Fields from 
Electron Density (EMFED), generating a continuous vector potential A (r, t) from phenomenological 
relations combining the unperturbed electron density of the material system, the material constants 
and the laws of optics. As an application of this model, we calculate in a time dependent approach the 
transition probability and the induced current density between the last bulk state below the Fermi 
energy and the first image state of a Cu(001) metallic surface. These observables are significantly 
modified by the spatial variation of the vector potential at the surface. The Coulomb gauge condition 
(V • A = 0), fullfilled everywhere else, breaks down near the surface. The difference between the s- 
and p- polarizations of the laser field partially unravels this effect. 

PACS numbers: 78.47,+p, 78.68.+m, 78.40.Kc 
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I. INTRODUCTION 

Because of its importance for the photoelectric effect, 
many studiesi^i&i have modeled the spatial variation of 
a laser electric field at a gas-solid (metal) interface in the 
long wavelength (LWL) domain (wavelengths A ^ 100 A 
or, equivalently, energies Hui ^ 124 eV). This electric 
field has been used to calculate the laser-matter inter- 
action and the associated photoelectric probability. In 
1975 Fcibclmanii has evaluated the spatial dependence of 
the laser vector potential for a gas-solid Jcllium interface 
solving the Maxwell equations using a nonlocal conduc- 
tivity tensor. This model has been applied by Levin- 
son et al& to photocmission and the results compared to 
experimental measurements. 

More recently, Miller et al£ have used the one step 
model of Mahan^, built on a Jellium potential displaced 
towards the vacuum to take into account a surface state, 
to study the direct, indirect and surface one photon pho- 
tocmission of the Ag(lll) surface. As in the work pre- 
sented in this paper, this model uses a laser-matter oper- 
ator in the velocity gauge written as a sum of a standard 
A ■ V and a surface (V • A) interaction term. Instead 
of an explicit calculation with a vector potential depen- 
dent on the spatial coordinate, the surface term has been 
fitted to the experimental spectrum to explain the un- 
symmetrical line shapes present in these data. Other 
models, derived from the work of Miller et al., have used 
the velocity gauge in the analysis of a two-photon pho- 
toemission experiment of metallic silver— near the surface 
plasmon resonance. 

Using the same velocity form of the laser-matter inter- 
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action, Tergiman et al& have calculated the contribution 
to a photo-current of a multiphoton excitation by solving 
analytically the Schrodinger equation for the electrons. 
These authors have modified the Fermi-Dirac distribu- 
tion describing the initial state of the electrons to in- 
clude the contributions from electron-electron collisions. 
With this model, Tergiman et al. were able to evaluate 
the linear and nonlinear multiphoton contributions to the 
photo-current as a function of the laser fluence. 

In recent years many experimental studies in photo- 
electron emission have used the versatile two-photon pho- 
tocmission (TPPE) techniquaiSiiiiiSii^ to study the tem- 
poral evolution of the electron dynamics and of the asso- 
ciated observables on clean and covered surfaces. Today 
the majority of theoretical simulations of TPPE are per- 
formed using the density matrix theory in its Liouvillc- 
von Neumann f orm 2ii2ii£ii§i±Li£. The density matrix for- 
mulation has the advantage of permitting the inclusion of 
the elastic and inelastic electron-electron collisions which 
constitute efficient desexitation channels when the phys- 
ical system is in interaction with the bath of the solid. 
Note that, except for the work of Pontius et al&, these 
density matrix formulations make use of the length gauge 
and the laser electric field is assumed to be independent 
of the spatial coordinate. 

In this article we discuss the laser-matter interaction at 
a gas-solid interface in the long wavelength domain. The 
sharp rise of the electron density at the gas-solid inter- 
face, taking place on a sub-nanometric scale, affects the 
electromagnetic fields through a steep variation of the re- 
fractive index even if the particular atoms of the solid are 
not directly "seen" by the electromagnetic fields. We try 
to answer the following question: what is the influence of 
the spatial dependence of the laser fields on the electron 
excitation at the interface and on the associated observ- 
ables? This question may be of importance especially 
in the presence of adsorbates and nano structures. In a 
formulation where the vector potential A is a function of 
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r, there is an additional contribution to the laser-matter 
interaction coming from the operator (V • A) mentioned 
above. This contribution is known as a "surface photo- 
electric effect" (see e.g. Desjonqueres and SpanjaarcUS). 

The above question is answered here in the framework 
of the Schrodingcr equation, for a clean metallic surface 
in the context of a single excitation between a discrctized 
state of the band of the solid and an image state. We have 
developed a phenomenological model, called Electromag- 
netic Fields from Electron Density (EMFED), where the 
electromagnetic fields arc explicitly function of the coor- 
dinate z normal to the surface. A simplified version of 
this model, where the electron density used to generate 
the vector potential was calculated from a Jellium solid, 
has already been published^SiSiiS. Here, the laser fields 
are calculated using macroscopic material constants, con- 
ductivity a, dielectric function e and refractive index n, 
which depend on the chemical nature of the bulk material 
and on the laser wavelength. These material constants, 
taken from experimental measurements tabulated by Pa- 
liliS&Si, are considered to be local relative to the spatial 
coordinate z, i.e. e(z,z') ~ s(z). To explicitly obtain 
the z-dependence of the laser field, the dielectric func- 
tion e of the metal is related to the unperturbed electron 
density p e {z) of the metallic surface evaluated from the 
corresponding wave functions of the occupied states of 
the system. As in the Jellium calculations performed by 
Lang and KohnSi using DFT (see also Jenning et alm^), 
the calculated electron density is normalized with respect 
to its value in the bulk. The dielectric function e(z) is 
then connected to the vector potential A (f, t) through 
its wave vector k ph via the refractive index h. 

Using the EMFED vector potential A (r , t) , the Hamil- 
tonian describing the interaction of the electrons of a 
metallic surface with the laser field is constructed. In the 
actual implementation of the model, the motion parallel 
and perpendicular to the surface are decoupled. Conse- 
quently, the wave function is expanded in terms of prod- 
ucts of parallel and perpendicular to the surface basis 
functions. In the directions parallel to the surface we 
use a linear combination of atomic orbitals which fulfills 
the Bloch periodicity condition, and the set of functions 
perpendicular to the surface is, for simplicity, restricted 
to a single resonant contribution. The time dependent 
Schrodinger equation is projected on this basis, yielding 
a system of coupled first order ordinary differential equa- 
tions (ODE) for the time variable. This ODE system of 
coupled equations is propagated for the duration of the 
laser pulse and the excitation probability and electron 
current density are calculated. 

Specifically we have studied a clean Cu(OOl) surface 
where the bulk continuum near the Fermi level is dis- 
cretized and we have considered a resonant excitation 
between this bulk Fermi and the first image states. The 
laser fluence used in the present simulation is in the range 
of the experiments by Velic et alSL for C$A§ on Cu(lll) 
or Kirchmann et alml for CqFq on Cu(lll), and is well 
within the perturbation regime. 



II. THE LASER FIELD AT THE INTERFACE: 
THE EMFED MODEL 

A. The vector potential of the laser field derived 
from the laws of optics 

The geometry considered here is depicted in Figure ^ 
the plane of incidence (POI) of the laser beam on the 
surface is defined as (xOz) 1 z being the direction normal 
to the surface, pointing towards the metal. The three- 
dimensional space is therefore divided in two parts: z < 
corresponds to the gas phase and z > corresponds to 
the solid. The angles Oi and 9t stand for the incident 
and transmitted angles. The angle between the POI and 
the crystallographic direction u of the surface is (p. The 
incoming electromagnetic wave, whose wave vector is de- 
noted by k ph , is p- polarized when the vector potential 
A (r, t) (notation P A) is in the POI, and s- polarized when 
this vector (notation S A) is normal to the POI. 
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FIG. 1: (color online) Schematic picture of reflexion and re- 
fraction of a polarized laser beam incident on a metallic (001) 
surface of a face centered cubic (FCC) crystal. The FCC lat- 
tice parameter do is shown and the shortest distance between 
two atoms on the surface corresponds to <2q = 0% = ao in 
the u and v directions. The vector potentials *A and P A, re- 
spectively for s- and p- polarizations, point in the directions 
perpendicular and parallel to the POI or xOz (see text for 
details) . 

In the following, we label (i), (r) and (t) the incident, 
reflected and transmitted components of the vector po- 
tential, taken as classical here. Taking into account the 
continuity of the fields at the interface (see e.g. Jack- 
son^), one writes the normal projection (coordinate z) 
of the vector potential of a p- polarized monochromatic 
laser beam as 



p A z (x,z,t) = sm6 t [M« + MM] (z < z P ) 
p A z (x,z,t) = sin0 t e{w,z) (z>z P ), 



(1) 



and the tangent projection (transverse coordinate x) as 
p A x (x,z,t) = COS0; IMW - p A^] (z < Zp) 



*A x (x,z,t) 



cos (9, [M« 
cos 8 t V A^ 



(z > z P ), 



(2) 



where zp is the position of the image plane which can be 
different from the position z = of the geometrical sur- 
face (see next subsection). In Eq. |JIJ, one should notice 



3 



the appearance of the complex relative local dielectric 
function £(oj, z), function of the z coordinate. 

For a 5- polarized monochromatic laser beam, the vec- 
tor potential is directed perpendicular to the POI, and 
therefore is parallel to the y direction in our coordinate 
system, and 

s A y (x,z,t) = *A« + *A^ (z < z P ) 
s A y (x,z,t) = *A® (z > z P ). 



(3) 



The incident, reflected and transmitted components of 
the vector potential used above are simply given by the 
standard classical expressions 

»°W(x,z,t) = A ei k * hx + kP * hz - ut ) 

* ol AW[x,z,t) = * ot R{u,z) Aoe^*-^*-"*) (4) 
»°W){x,z,i) = ?° l T(uj,z)Aoe< k * hx + k * hz -" t ) 

where pol corresponds to the p or s linear polarizations. 
The reflection pal R(cu,z) and transmission pal T(uj,z) co- 
efficients are given by 



and 



p R(oj, z) 
v T{u:,z) 

s R(lu, z) 
*T(u,z) 
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(5a) 
(5b) 

(6a) 
(6b) 



cos Oi + n(u), z) cos t ' 

where the complex refraction index h(u>, z), also function 
of z coordinate, is obtained from the EMFED model as 
described in the next subsection. The attenuation of the 
field in the metal is simply taken into account through a 
complex wave vector 

F^Su/c. (7) 

In Eq. (@J, the normalization factor Aq defines the am- 
plitude of the incident laser vector potential. It can be 
calculated from the fluence F or the intensity 3q of the 
laser pulse. If the vector potential is described by a sin 2 
pulse shape of full width at half maximum (FWHM) r, 
one gets 



A)(a.u.) = 



1.66 10 



-3 



'F(J/m 2 ) 



(8) 



i w(a.u.) y i~(fs) 

The numerical relation between the fluence F and the 
laser intensity 3q for this pulse shape is 

F(J/m 2 ) = 0.9708 3 (W/m 2 ) r(s). 

The angle of refraction t is, in this phenomenological 
model, a complex function of the incident laser angular 
frequency to and of the coordinate z. It is calculated from 
the law of Sncll-Descartcs 



sin 



sm ft t = — 



h(u), z) 



(9) 



B. Material constants from the conductivity and 
the electron density 

The equations llll-'>l) are standard expressions of conti- 
nuity at an abrupt interface. Here the material constants, 
as for example the dielectric function e or the refractive 
index h, are functions of the coordinate z. This depen- 
dence also induces a dependence relative to z of the as- 
sociated wave number k ph (Eq. J7J)) and of the refraction 
angle t (Eq. ©). 

Following Drude's model 2 ^, at optical and higher fre- 
quencies, both bound and conduction electrons of a metal 
contribute to the dielectric function through the relation 



:(lo) = e b (uj) + i 



(10) 



where £b(u) is the relative permittivity due to the bound 
electrons, £o denotes the vacuum permittivity (only non- 
magnetic materials are considered here with a relative 
permeability \i — 1) and &(u>) is the complex conductiv- 
ity (see e.g. Jackson 2 ^, chapter 7). The Drude model 
also gives the following analytical expression for the con- 
ductivity 



m* (70 — i lu) 



(11) 



where e is the charge of the electron, to* its effective 
mass, p s the electron density of the bulk and 70 the elec- 
tron friction coefficient. 

The above equations relate e(uj) and p s and can be 
used as a template to model in a simple way the spatial 
dependence of the dielectric function at the surface 



z) = 1+ (i s {u) - 1) p e {z), 



(12) 



where e s (tu) is the complex relative dielectric function in 
the solid and p e {z) is the relative electron density normal- 
ized with respect to the bulk density p s . Equation itT^) 
has the advantage of imposing a continuous variation of 
e(w, z) through the interface with a linear dependence 
on p e (z). In addition this expression shows the correct 
asymptotic limits 



-00 j 
-00) 



1 

e s (w) 



Equation (fT2^l is the heart of the EMFED model which 
uses both the unperturbed relative electron density of 
the material system p e (z) under study and the tabulated 
values^Si of the associated complex refractive index of 
the bulk (n 2 (w) = s s (uj)). 

Using the density functional theory (DFT), the spatial 
variation of the electron density p e [z) at a gas-Jellium 
interface has been calculated by Lang and Kohn^ for 
several electrons densities in the bulk p s , or Wigner-Seitz 
radii r s = (3/(47rp s )) 1 / 3 . A calculation of this electron 
density from the occupied states of the material system, 
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using a smooth model Jcllium potential at the inter- 
face due to Jennings et alm&, gives similar results. In 
the present implementation of the EMFED model to a 
metallic surface, we calculate this electron density from 
a more accurate model potential perpendicular to the 
surface parametrized by Chulkov et al£&. The Jellium 
and Chulkov et al. potentials are displayed as red solid 
and black dashed lines in the lower part of Figure [21 F° r 
a given potential, the total relative electron density p e (z) 
is expressed as a discretized sum 

Pe {z) cc Y, I^WI 2 ( 13 ) 

I—l, np 

running over all the occupied states up to the Fermi level. 
r)i(z) denote here the wavefunctions associated with 
the eigenstates of the one-dimensional time-independent 
Schrodinger equation in the z direction. This relative 
electron density p e (z) is normalized with respect to the 
mean average value of the density p s deep in the bulk. 
The sum is discretized because we have adopted a dis- 
crete variable representation (DVR) approachi^S for the 
description of the continuum band structure of the metal. 
The electron densities calculated with the potential of 
Chulkov et al. and with the Jellium potential are dis- 
played as red solid and black dashed lines in the upper 
part of Figure |3 For the Cu(OOl) surface studied here, 
the electrons at the Fermi level belong to the occupied 
band, and, because of the use of the DVR, the band is 
discretized. The last discrete state below the Fermi en- 
ergy, hereafter called bulk Fermi state, and the first image 
state (n = 1) are displayed in the lower part of Figure |21 
These states participate in the photoexcitation consid- 
ered here. Note that in the present case these states and 
the total electron density extend over several angstroms 
in the vacuum. Note also the difference in the amplitude 
of the oscillations of the densities and potentials between 
the Chulkov et al. and Jellium models. This difference 
is due to the explicit inclusion of the atomic structure in 
the first but not in the second model. 

The refraction index is obtained from the standand 
relation to the dielectric function 

n(w, z) = y/e(w, z) . (14) 

Using the electron density (|13|l . the analytic form of the 
dielectric function l|12f) . and the associated refraction in- 
dex one can calculate, from the expressions of subsec- 
tion [^2J the vector potential A (r, t) of the laser field at 
the interface. This simple procedure allows us to intro- 
duce the spatial variation of the vector potential at the 
gas-solid interface. 

The EMFED model explained above can tentatively be 
understood from the point of view of the solutions of the 
Maxwell equations. Far in the vacuum the asymptotic 
condition corresponds to the incident plane wave p 
of Eq. J3J. The propagation of the Maxwell equations 
towards the interface smoothly accumulates the presence 
of the electron density in the vector potential. At the 



end of the propagation far in the bulk the vector poten- 
tial corresponds to the transmitted component f^A^. In 
the EMFED model, the reflected component contributes 
near the surface, therefore allowing for the continuity of 
the electromagnetic fields at the interface. This result 
is therefore qualitatively similar to a model of Feibel- 
manS who solved numerically Maxwell equations for a 
vacuum/ Jcllium system. 
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FIG. 2: (color online) Lower graph (a): The electron-matter 
interaction potentials for Cu(001) using the Chulkov et al. 
(red solid line) or Jellium with r 3 —2.67 a.u. (black dashed 
line) models. The probability densities |t7j(.z)| 2 associated 
with the bulk Fermi and first image states accommodated 
by the Chulkov potential are also shown at their energy (- 
0.12 eV from Ef; -0.574 eV from the vacuum) in blue and 
green respectively. Upper graph (b): the total electron den- 
sities calculated using Eq. I I -'it with the same colors and line 
types. The vertical line at zp = —1.31 A corresponds to the 
position of the image plane P at p e (zp) — 0.5 where the con- 
tinuity of the vector potential is enforced. 

For the s- polarization, the vector potential presents 
a single y component displayed in Figure [3J It shows a 
smooth variation everywhere and particularly near the 
image plane of the Cu(001) system at zp ~ —1.31 A cor- 
responding to the electron density p e {zp) = 0.5 where, 
following Eqs. Q), @ and ©, the continuity of the vec- 
tor potential is enforced. This image plane can be shifted 
from p e {zp) = 0.2 to p e (zp) = 0.8 with nearly no con- 
sequence on the vector potential itself. Similar smooth 
behavior is obtained in the calculation of Fcibclmar. 33 . 
The dashed lines in Figure [21 correspond to a calculation 
with a constant vector potential in each medium where 
only the incident and transmitted components are taken 
into account. 

Various derivations and models have also used the 
material constants in the two media together with the 
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continuity at the interface to obtain the corresponding 
properties of the laser-matter interaction. For example, 
Feibelman24 obtained in this way the surface plasmon 
frequency and Liebsch and Schaiehs* explained the silver 
anomalous dispersion of the surface plasmon frequency. 

Because we are not solving explicitly the Maxwell 
equations for the electromagnetic fields, the problem of 
the gauge used here has not been touched up to now. 
How behaves the external electromagnetic field in the two 
media (vacuum and bulk) and at the interface from the 
point of view of the gauge of the electromagnetic wave? 
The vector potential can be decomposed in perpendicu- 
lar and parallel components with respect to the surface 
plane as A = A± + An. A simple look at Figure Ogives 



p A, = 



P A 



for the p- polarization, and 



'A ± 
S A« 



p A z (x,z) 
p A x (x,z) 





S Ay (X, Z) 



(15) 



(16) 



for the s- polarization. In this last case of s- polariza- 
tion, it is already clear that V ■ S A — everywhere, and 
the Coulomb gauge condition is therefore verified in our 
model for this specific polarization. The case of p- po- 
larization is clearly different since the fast variation of 
Mj_ with z in the vicinity of the surface implies that 
V • P A ^ when approaching the metal. As a conse- 
quence, the electromagnetic fields given by the EMFED 
model do not verify the Coulomb gauge for the p- polar- 
ization. This discussion also applies to the Schrodingcr 
equation, and the laser-matter interaction term in the 
Hamiltonian includes a non-Coulomb gauge term pro- 
portional to (V • A). 



III. A TIME DEPENDENT MODEL FOR THE 
ELECTRONIC DYNAMICS AT THE INTERFACE 

A. The Hamiltonian 

The total electronic Hamiltonian including the laser- 
matter interaction is similar to the Hamiltonian for the 
unperturbed system except that the electron momentum 
p is replaced by (j> + gA), where e and A(f,t) are the 
charge of the electron and the vector potential of the 
laser field. The total Hamiltonian reads 

H{r,t) = H a (r)+H L (f,t)+H Q (r,t), (17) 

where the zero order Hamiltonian is 

r2 



(18) 
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FIG. 3: (color online) Real and imaginary components y of 
the vector potential A as a function of the coordinate z for a 
s- polarized laser beam. Following Eq. the x dependence 
in Eq. (|160 corresponds to a plane wave and therefore the 
dependence of the vector potential relative to this variable is 
not shown in the present figure. The photon energy and the 
laser fluence are 4.169 eV and 10 /^J/cm 2 respectively. The 
real and imaginary refraction indexes of the bulk metallic Cu, 
n=1.386 and /t=1.687 at 4.2 eV, are taken from the tables of 
Palik22*24. The black Re[A] and red Im[A] solid lines are the 
results obtained with the EMFED model, while the dashed 
lines correspond to the incident and transmitted fields in a 
simple approach from optics with constant vector potentials 
(quotted as gas/solid in the figures below). 



can be decomposed as 



he 



H L (r,t) = — [2A.V + (V.i)] 
H Q (r,t) = ^[A-A 



(19a) 
(19b) 



where the parentheses in the first equation above mean 
that the operator V acts on the vector potential A(r, t) 
only. The superscripts L and Q stand respectively for 
the linear and quadratic terms in A. 



TABLE I: Decomposition of the laser-matter interaction 
Hamiltonian for a p- and s- polarized vector potential, re- 
spectively p A(f, t) and s A(f,t). 
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Ay 





The laser-matter interaction terms of the Hamiltonian 



The analysis given in the preceding section has re- 
vealled a non-zero divergence of the vector potential 
(V • A) in the case of p- polarization. This non-Coulomb 
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gauge contribution develops in two terms, one perpen- 
dicular (dA z /dz), and the other parallel (dA x /dx) to 
the surface (see Table QJ. These two terms, called sur- 
face terms, give rise to the so called surface photoelectric 
effect (see e.g. Desjonqueres and SpanjaarcUS). For s- 
polarization, these two terms do not appear, but Table [I] 
shows that two standard contributions remain in the in- 
teraction Hamiltonian: 2 A y d/dy and Ay. 

Table H] also shows that a characterization of the laser 
field can be performed by calculating the difference be- 
tween the observables measured for p and s- polariza- 
tions. For equivalent directions x and y parallel to 
the surface, the contributions from A x + 2 A x d/dx of 
p and A y + 2 A y d/dy of s- polarizations are equiva- 
lent and eliminated from the interaction. Unfortunately 
A\ +2 A z d/dz is not eliminated and therefore one can- 
not measure directly the influence of the surface terms 
(V • A) = {dA x /dx) + (dA z /dz). 



B. The time dependent Schrodinger equation 

The time-dependent Schrodinger equation describing 
the interaction of the laser field with the electrons at the 
gas-solid interface reads 

ih%-9{?,t)=H(?,t)*(?,t), (20) 
at 

where ^(r, t) denotes the electronic time-dependent 
wavefunction. This wave function is expanded on a series 
of time independent basis functions <&/ (r ) 



*(f,i) = 5^$j(f)(7 z (*) = #C, 



(21) 



where Ci(t) are the time dependent expansion coeffi- 
cients. The second expression corresponds to a matrix 
notation. 

Figure shows that the presence of the two media 
breaks the periodicity in the direction z normal to the 
surface. On the other hand, and since there is no elec- 
tron density in the vacuum, we will treat the system as a 
periodic system in the parallel direction even in the vac- 
cum. In the present version of the model we write the 
basis wave function &i(r) of Eq. (|21|l as a product 



$/(r) = r)i(z) //(f|i 



(22) 



where i]i(z), already defined with Eq. (|13f) . is a 
wavefunction associated with the eigenstates of the 
one-dimensional time-independent Schrodinger equation 
along the perpendicular direction z. Along the surface 
plane direction r\\ = (a;, y), the basis functions fi{r\\) are 
written as 



//(ni) = E 



Xl{f\\ -Rl): 



(23) 



where denotes the parallel electron wave vector. In 
this expression, xi denotes a function defined in the el- 
ementary cell. The vectors Rl allow for the translation 
of this localized wavefunction in the parallel direction. 
This generates a two-dimensional periodic function in the 
surface plane which fulfills the Bloch periodicity condi- 
tion (see for instance Ziman^i) . This approach is similar 
to the method of linear combinations of atomic orbitals 
(LCAO) for a solid. 

Inserting the expansion (|21() in the time-dependent 
Schroringcr equation H2U|) . and projecting it on the time 
independent basis functions $/(r) yields a coupled sys- 
tem of first order ordinary differential equations (ODE) 
for the expansion coefficients Cj{t) 

iK^-=Y,{$i\H{r f ,mj)Cj{t). (24) 
j 



dt 



The diagonal terms can be written as 

($i\H{?, t)\$i) = E± +E\ + {^!\H l + , (25) 

where Ej- and Ej denote the energies associated with 
the wavefunctions r)i(z) and fi(r\\) respectively. 

For low laser intensities, the off-diagonal laser-matter 
interaction matrix element (&i\H L + H®\$j) can be ap- 
proximated by its linear component Here 
this term is an integral over the entire space. On the 
other hand, the current density, Eq. I|30b(l below, also 
contains this integral but taken over a restricted domain 
in the perpendicular direction, namely from — oo to z. 
One can write this integral, to be used in these two ex- 
pressions, as 



he 
2im 



Ju{z) e 



(26) 



where the calculation of Ju(z) is detailed in Ap- 
pendix [S] 

Using the model explained above, two observables are 
calculated: the excitation probability and the electron 
current density. The population of the different eigen- 
states and therefore the excitation probability are ob- 
tained from the overlap between the states $ / (r ) and the 
final wave packet ^>(f,tf) at the end of the pulse (t = tj) 



P I (t f ) = \($ I (r) \ *{r,t f ))\ 2 = \C I {t f )\ 2 , 



(27) 



The expression for the electron current density, includ- 
ing the linear momentum of the electron p and of the laser 
field eA, can be written using the kinetic energy operator 
T = (p + eA) 2 /2 m (see e.g. Messiah^) as 



in 



dz' I dS 

'See 



* f f * - (f^) t ^ , (28) 



where S ce u denotes the surface of the elementary cell 
along the parallel direction. This total current can be 
split in two parts 



I(z,t) = I (z,t) + I int {z,t). 



(29) 
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with 



lim 



I dz'l dS (*tV 2 * - V 2 ^*) 



and 



2m 



-ct 



' •fint(z) 



C (30a) 



(30b) 



As in Eq. (|2ip. in the previous equations $ is the line 
vector of the basis functions and C is the column vector 
of the coefficients. The term 3i nt (z) is a matrix whose 
elements are the integrals Jjj{z) given in Eqs. l|A4a(l and 
(|A4b|) . In the calculation of Tint , the quadratic terms in 
A have been neglected. 



IV. THE PHYSICAL MODEL AND THE 
EXCITATION DYNAMICS 

A. Details of the physical model 

In the present model calculation we concentrate on the 
influence of the spatial dependence of the laser field in the 
long wavelength domain fsection fll Afl on the observables 
associated with a photoexcitation process. 

The model representing the material system is re- 
stricted to its simplest mcaningfull form which is sum- 
marized below. First, the clean Cu(001) surface has been 
selected since it corresponds to a simple surface structure. 
The vectors u and v, used to generate the (001) surface 
displayed in Figure^ are not the primitive vectors of the 
FCC solid and the associated potential V(u,v,z) is not 
separable. For example, in the plane of incidence (xOz) 
there is a shift by ao/2 in the z direction between two 
adjacent columns of atoms (see the left part of Figure^). 
Starting from a DFT calculation, Chulkov et al^SL have 
obtained a model potential in the z direction where the 
periodicity is ao/2. With this potential in the z direction, 
the tridimensional potential V(u,v,z), now separable in 
the first approximation, is used in the present simulation. 

Secondly, solving the time-independent Schrodinger 
equation in this potential yields the set of wave func- 
tions defined in Eqs. 1221) and (|23fl . If one discrctizes 
the bands of the solid using a DVR approach^iiS, the 
Chulkov potential for the e-Cu(OOl) system accommo- 
dates discrctized bulk and image states. In the z direc- 
tion our grid extends from 70 A in the vacuum to the 
200 A in the bulk, with 5881 grid points. As the initial 
state we select the last occupied discretized bulk state be- 
low the Fermi energy located at -0.12 eV from the Fermi 
energy already designated under the name of bulk Fermi 
state. The final state is the n = 1 image state located at 
-0.574 eV from the vacuum energy. These two states are 
separated by 4.169 eV and their densities are presented 
in the lower part of figure [21 

In the direction parallel to the surface we use a peri- 
odic wave function given in Eq. (12.' j|) . We do not perform 



an explicit calculation to obtain the optimized wave func- 
tions. Instead, in the elementary cell, the localized wave 
functions xi{ u i v ) are two-dimensional harmonic oscilla- 
tor wavefunctions. Following the usual selection rules, we 
have restricted this basis set to the s and p symmetries 
only 



Xs{u,v) = 




(31a) 


X P „(u,v) = 




(31b) 


X Pv (u,v) = 




(31c) 



As shown in Eq. (|23f) . the periodic wave function fi(r\\) 
is obtained from these localized basis functions translat- 
ing them by a multiple of and Oq (see the caption 
of Figure G). On the Cu(001) surface this distance is 
a o = a W = 2.556 A and the symmetry of the surface cut 
imposes the use of the same exponent £ in the u and v 
directions. This exponent was chosen such that the re- 
sulting wave functions are confined in the elementary cell. 
For simplicity, and because no optimization is performed 

parallel to the surface, the energies Ej of the s, p u and 
p v states in Eq. (|25J) arc set to zero. 

The vector potential A(r,t) has been obtained using 
the EMFED model following the prescriptions given in 
section [H] where the electron density is calculated us- 
ing the Eq. (|13fl and the wave functions rji(z) (Eq. i|22|) 
are obtained using the DVR calculations performed using 
the grid given above. At the photon excitation energy of 
4.169 eV, the complex refraction index h = 1.386+1.687 i 
of the solid is taken from the tables of Palik 2 ^. It is used 
to calculate the dielectric constant i s (ui) of the solid and, 
combined with the electron density from Eq. (|13fl . the 
dielectric function e(u>, z) in Eq. (|12[> . Finally, the refrac- 
tion index n(u>, z) is calculated from Eq. (|14|) and used 
to construct the EFMED vector potential in subsection 

nra 

In the calculations, the incident angle of the laser beam 
is fixed at 6*^=45 degrees except for the one presented in 
figure |21 where this angle is varied. For the Cu(001) sur- 
face, the observables calculated in the present work are 
independent of the orientation angle (p of the POI (see 
Fig. GJ. We have therefore fixed the orientation of the 
plane of incidence (POI) parallel to the u axis of the ma- 
terial system i.e. tp=0. With these parameters, the vector 
potential A(r,t) is calculated using Eqs. (lll.'ill . We con- 
sider a pulsed laser with a peak intensity of 10 s W/cm 2 
and a temporal FWHM of r = 80 fs resulting in a flu- 
ence of 10 /xJ/cm 2 . This fluence is well within the per- 
turbation regime, where space charges near the electron 
analyzer do not influence TPPE measurementsiSiSi. 

The calculations are performed at the T point and the 
dispersion is neglected. Therefore, in equation (|A12I) of 
the appendix, only the first term is calculated. With the 
present choice of basis functions parallel to the surface, 



8 



the laser-matter interaction terms between the s and p u 
or p v states, given in Eq. (|A13() . reduce to 



B. Excitation dynamics 



UT) 00,00 

U IJ 



d 

X * s du Xpu du dv = ^ ' ^ 32 ^ 



where s and p u belong to different states / perpendic- 
ular to the surface. The ODE system of coupled equa- 
tions l|24[) contains six states, three s, p u and p v for the 
initial and the final states. In the unambiguous cases 
these initial and final states will be referred simply to 
as i and /. The integrals over the basis functions for 
the e-surface potential and the laser-matter interaction 
have been obtained once for all and the system of six 
coupled ODE (|24(l has been solved at each time step by 
a computer program based on a predictor-corrector algo- 
rithmffi. The obtained time dependent expansion coeffi- 
cients permit the calculation of the excitation probability 
(Eq. J2Z}) and of the current density (Eqs. l^ - pp l. 
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FIG. 4: (color online) Temporal evolution of the excitation 
probability at resonant energy huj — Ef — Ei — 4.169 eV for a 
pulsed laser FWHM of r = 80 fs and intensity of 10 s W/cm 2 . 
The displayed time interval corresponds to a total duration of 
220 fs. Upper graph: Excitation probability calculated with 
the EMFED vector potential with the laser polarizations p 
and s. The direction of the incident light is Si =45 degrees and 
the POI is parallel to the (100) direction that corresponds to 
tp = in Figure Lower graph: Excitation probability calcu- 
lated with the p-polarization for the EMFED vector potential 
(black full line), for an abrupt vector potential at the inter- 
face (gas/solid, blue dashed line) and for a constant vector 
potential from the gas phase (gas, green dotted dashed line). 




FIG. 5: (color online) Current density 7o, corresponding to 
the kinetic energy operator V 2 obtained from the Eq. 130at . 
as a function of the coordinate z at the peak of the laser pulse. 
This current density is normalized to the maximum of s Iq us- 
ing the normalization factor of 6.5 xlO . The characteristics 
of the laser are given in the caption of Figure 2] 
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FIG. 6: Interaction current density Ii nt for the p polarization 
of the laser corresponding to the operator 2 A ■ V + (V • A) 
at the peak of the laser pulse as a function of the coordinate 
z calculated from the Eq. I30bl . This current density is 
normalized using the factor 3.37xl0 -5 . The characteristics 
of the laser pulse are given in the caption of Figure [1] 



The excitation probability between the bulk Fermi 
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FIG. 7: (color online) Temporal evolution of the extrema of 
the current density Io at resonance hui — 4.169 eV for the p- 
and s- polarizations of the laser. The characteristics of the 
laser pulse are given in the caption of Figure 2] The cur- 
rent density is normalized to the maximum of s 7o using the 
normalization factor 1.31 xlO -4 . Upper graph: the displayed 
temporal region corresponds to the peak of the laser pulse. 
Lower graph: envelope of the temporal evolution of the ex- 
trema of the current density To- 



state and the first image state has been calculated at 
resonance for the two laser polarizations p and s. The 
evolution of this excitation probability during the laser 
pulse is presented in Figure 0] It is the sum of the ex- 
citation probabilities of the s, p u and p v image states 
calculated using the Eq. (|27[) . This observable does not 
display the oscillatory pattern related to the phase of 
the laser field since the square in eq. (|27[l annihilate the 
phase of the expansion coefficients and the rise and the 
stabilization of the excitation probability is smooth. 

As seen in the upper part of Figure 0] the excitation 
probabilities of the s and p- polarizations are different. 
The expression of the excitation probability for the s po- 
larization contains a single contribution from the A y com- 
ponent of the vector potential. For p polarization both 
A x and A z contribute, giving rise to a destructive inter- 
ference which explains the present lower probability for 
this polarization. 

In the lower part of Figure 0] one compares the exci- 
tation probability for the p polarization calculated with 
the vector potential of the EMFED model, with the one 



calculated with a discontinuous vector potential at the 
interface (label gas/solid, see also Fig.^J. A comparison 
is also performed with a vector potential independent of 
the spatial coordinate (label gas). The excitation proba- 
bility calculated using the gas/solid discontinuous poten- 
tial gives a result 25 % lower than that calculated with 
EMFED. The simulation using the spatially independent 
constant vector potential gives an excitation probability 
24 times smaller than that obtained using the EMFED 
model. Taking into account the spatial variation of the 
laser field at the interface is therefore crucial for a realis- 
tic description of the electron dynamics at the interface. 

Moreover we have calculated the current density de- 
composed into two contributions: Iq and / m t. The cur- 
rent density Io is three orders of magnitude larger than 
lint and therefore Iq presents all the characteristics of 
the total current density. The spatial dependence of the 
current density Iq and Ii„t at the peak of the laser pulse, 
obtained from the expressions (|3Ua[) and (|30b(l ). is pre- 
sented respectively in Fig. and 

The Iq current density displays spatial oscillations in 
the z direction that mimic the oscillations in the electron- 
solid interaction potential shown in the lower part of 
Fig. |21 The Ii n t current density shows similar but atten- 
uated oscillations. The Eqs. (|3Ua|) and l|30b(l show that 
the Iq and Ii n t current densities are calculated as defi- 
nite integrals from — oo to z. Changing the upper limit 
z adds positive or negative contributions to the current 
giving rise to the observed oscillations. As for the exci- 
tation probability in Figure 01 and for the same reasons 
the s current density is larger than the p current den- 
sity. Outside the surface region the current density is 
zero therefore at a nanoscopic scale there is no charge 
migration. 

The lower part of Fig. [31 displays the Iq current density 
for the different approximations of the vector potential. 
One sees a similar behavior to the one obtained for the 
excitation probability. These two observables are of dif- 
ferent nature and one concludes that the introduction of 
the variation in space of the vector potential is essential 
for a proper modeling at the gas-solid interface. 

Figure presents the current density Z;„ t for the p po- 
larization of the laser. For the s polarization the corre- 
sponding current density Ii nt is zero and is not displayed. 
This result can be explained as follows. For the p- po- 
larization, the numerical derivatives of the basis func- 
tions, relative to the z coordinate, do not compensate in 
Eq. (|30bf) resulting in a non-zero interaction current den- 
sity Ii n t ■ The zero contribution for the s- polarization is 
traced back to the integral Eq. (|30b|l where in this case 
the two terms compensate each other due to the con- 
servation of the momentum parallel to the surface (see 
appendix) . 

Next we present in Figure [71 the temporal evolution 
of the extrema of Iq, around z ~ -2 A. The upper graph 
presents a small region of the temporal evolution of these 
extrema near the peak of the pulse for the p- and s- po- 
larizations. As expected, the current density oscillates 
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FIG. 8: (color online) Energy behavior near the resonance en- 
ergy at huj = 4.169 eV of the excitation probability for t — > oo 
(a) and of the maximum of the current density To (b) approxi- 
mately located in the vacuum at z ~-2 A. This current density 
is normalized to the maximum of p 7o using the normalization 
factor 6.8xl0 -5 . The characteristics of the laser pulse are 
given in the caption of Fig. 0] 



at the Bohr frequency lob = (Ef — Ei)/K ~ 6.3 fs _1 . 
The calculations for the different polarizations appear 
in the same order as the excitation probability. More- 
over, there is a small temporal dephasing between the 
p- and s- polarizations. A possible origin for this time 
shift is the destructive interference between the P A Z and 
V A X contributions to the current (see Eqs. J2HJ), l|A4a|) 
and \ A4b|) ) . In opposition to the excitation probability 
where the laser phase is lost, the observed oscillatory be- 
havior in Iq is due to the relative phase of the products 
of the complex expansion coefficients C* Cf or C'j Ci in 
Eq. (|30a|) . The envelope of the evolution of these extrema 
of the current density Iq , presented in the lower part of 
Figure displays a smooth variation similar to the evo- 
lution of the excitation probability (Figure QJ. Because 
the interaction current Zj n j is about three order of magni- 
tude smaller than Iq (see Fig. GD and EI), we do not present 
the similar temporal variation of this current. 

A final series of results is shown in Figures |H1 and [5] 
The first figure displays the behavior of the extrema of 
the calculated observables as a function of the photon 
energy near the resonant excitation energy between the 
bulk Fermi state and the n = 1 image state. The excita- 
tion probability at t— > oo and the maximum of the cur- 
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FIG. 9: Excitation probability for t — > oo (a) and maxi- 
mum of the current density To (b) located in the vacuum 
approximately at z ~-2 A and normalized with the factor 
1.2xl0 -4 as a function of the incident angle 8i of the laser 
beam. The calculations are performed at the excitation en- 
ergy hco — 4.179 eV and the characteristics of the laser pulse 
are given in the caption of Figure [1] The calculations cor- 
respond to the laser polarizations p (full black line) and s 
(dashed red line). 



rent density Iq display the usual behavior where a max- 
imum appears at the resonance energy. The FWHM of 
these Lorentzian-like curves for the excitation probability 
and for the current density are 0.026 eV and 0.035 eV, 
in agreement with the bandwith of the pulse. 

The second figure displays the evolution of the excita- 
tion probability and of the maximum current density Jo 
with the incident angle 0i of the laser beam (see Fig.^l. 
As one expects, when the incidence is normal to the sur- 
face, the observables have the same value for the p- and 
s- polarizations. At grazing incidence the laser field does 
not penetrate in the solid and no excitation occurs. But 
as soon as the incident angle Ot diminishes, for example 
for Bi = 89 deg these observables rise sharply. 



V. DISCUSSION AND CONCLUSION 

This paper presents a detailed analysis of the laser- 
matter interaction at the gas-solid interface in the long 
wavelength domain, i.e. A ^ 100 A. The main goal of the 
present work is the identification of the role played by 
the spatial variation of the laser fields, due to the sudden 
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rise of the electron density at the surface, on the electron 
dynamics at the gas-solid interface. 

First let us briefly recall the model. For a one photon 
excitation, our theoretical model calculates the excita- 
tion probability and the associated electron current den- 
sity between a state of the solid near the Fermi energy 
(called bulk Fermi state) and the n = 1 image state. The 
laser-matter interaction Hamiltonian (subsection IIII A|) 
is calculated using the vector potential of the laser field 
obtained from the EMFED model as a function of the z 
coordinate (section [Hj. The electron dynamics at the in- 
terface is obtained by the solution of the time dependent 
Schrodinger equation (subsection IIII Bjl . 

For the Cu(OOl) surface, the present model makes use 
of a separable potential parallel and perpendicular to the 
surface (see subsection IIV A(l . the last component of the 
potential being taken from the work of Chulkov et al~. 

From the total unperturbed electron density p e {r) of 
the material system, the EMFED model obtains, using 
the constants of the material, a realistic vector potential 
A(f, t) function of the spatial coordinate for an arbitrary 
incident angle of the photon. The only needed ingredient 
is an electron-matter interaction potential. 

The total wave function is expanded (Eq. I^lf) ) in a se- 
ries of six discrete basis functions (see subsections IIII Bl 
and IIV Afl . The first three are constructed as a simple 
product of the ground bulk Fermi state in the perpen- 
dicular direction and asorp orbital parallel to the sur- 
face. The last three are constructed similarly but with 
the n = 1 image state in the perpendicular direction. 
The projection on these basis functions of the time de- 
pendent Schrodinger equation gives rise to a system of 
coupled first order ordinary differential equations in time 
(Eq. I|24|) ') which are solved for each time step by a pre- 
dictor corrector algorithm. 

Next, the results of the present model can be sum- 
marized as follows. The analysis of Figures [3] and 
concludes that, near the gas-solid interface, a realistic 
model of the electron dynamics in the presence of the 
laser-matter interaction requires the inclusion of the spa- 
tial variation of the laser field even in the long wavelength 
domain. 

The excitation probability (Fig. 0} and the current 
density (Fig.|3J are different for the p and s- polarizations 
of the laser, giving an indication on the so called surface 
photoelectric effect appearing in the neighborhood, ap- 
proximately ± 3-5 A, of the interface. Precisely, Table [I] 
gives the different contributions to the interaction Hamil- 
tonian and one sees that, even if the vector potential 
projections x and y are equivalent, the x scalar product 
contains the derivative with respect to x of the vector 
potential whereas the y scalar product does not contain 
a similar surface term. Therefore, from the difference be- 
tween the observables for the p and s- polarizations, in 
the present excitation between discrete states or an exci- 
tation giving rise to an electron emission in the vacuum, 
one cannot estimate precisely the surface photoelectric 
effect. 



The current density Io(z,t) (Fig. [3J) reproduces the 
spatial oscillations of the electron-solid potential as well 
as the temporal oscillations of the laser field. Contrary 
to the excitation probability, this observable is unravel 
the phase of the time evolution of the laser. The current 
density Ii nt (Fig.[f)|) is directly related to the laser-matter 
interaction term in the Hamiltonian and, for photoemis- 
sion, to the surface photoelectric effect. But this term is 
weak and cannot be measured experimentally 

Finally, let us discuss the approximations and the 
limitations of the present model. First, the separable 
electron-surface potential contains a contribution in the z 
direction explicitly calculated and a contribution parallel 
to the surface that is not explicitly calculated. But this 
separable potential could be replaced by a more realis- 
tic tridimensional potential obtained using, for example, 
a DFT calculation. In addition, the present simulation 
of the excitation between discrete states can easily be 
extended to the photocmission and, because our wave 
function is expanded on a basis set, the number of these 
basis functions can easily be raised. 

Secondly, the use of tabulated material constants in 
our model implies that this model is able to account for 
the presence of the plasma oscillations in the bulk as one 
of us has shown for a Jellium Al(100) surface 22 . One 
can also account for the change in the observables at the 
surface plasmon frequency because in the present model 
the dielectric function depends on z. The dependence 
introduced here is similar to the one used by Feibelman^i 
when he derives the frequency of the surface plasmon. 

Let us now replace the present work in a more gen- 
eral context. The consequences of the spatial variation 
of the vector potential have been highlighted through a 
numerical calculation of the excitation probability and of 
the current density. But other physical observables qual- 
itatively sensitive to such a spatial variation of the vec- 
tor potential, have to be identified, calculated and mea- 
sured. For example, time resolved TPPE experiments, 
coupled with some interference measurements^ or combi- 
nation of TPPE and photocmission electron microscope 
(PEEM)iMi, could be efficient experimental techniques 
for observables sensitive to the spatial variation of the 
vector potential. Particularly the last cited technique 
TPPE/PEEM allows for an imaging of single nano ob- 
jects of a size of about 50 nm thus diminishing the aver- 
aging made over many objects in standard experiments. 
Combination of a laser with a STMM^ in a single simul- 
taneous experiment can be also used for measuring the 
observables related to these nano objects. 

One can easily extend the present model to single pho- 
ton photocmission but such a model remains oversimpli- 
fied for the previously cited experimental techniques. In 
fact, to model a TPPE experiment, one has to take into 
account, in addition to the direct excitation processes 
and associated continua of the electron, the electron- 
electron relaxation at the surface and in the bulk. Mod- 
cling such an experiment, implies the combination of the 
EMFED vector potential with a density matrix formula- 
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tion£ii2ii£ii£iiLi&. Our approach calculates some of the 
ingredients needed for such a density matrix model, as for 
example the unperturbed states and the interaction ma- 
trix elements with a spatially varying vector potential. 
Other needed ingredients, like the elastic and inelastic 
parameters of the electron-electron collisions correspond- 
ing to the relaxation to the bath, can be taken from the 
existing literature. 

Concerning other systems to be studied, different clean 
metallic surfaces can be selected, but adsorbate a 10 ' 12 or 
nano objects on surface s 11 ! 13 ! 41 could also be very sensi- 
tive to a spatially varying vector potential. The calcula- 
tion of the observables as a function of photon frequency, 
particularly at the bulk and surface plasmon resonance, 
is the next goal in our modeling. Generalization to semi- 
conductor surfaces is also of interest. However, for these 
surfaces reconstruction takes place and the wave func- 
tions parallel to the surface should be functions of the 
slab layer. 

In summary, the present model clearly shows that the 
spatial variation of the vector potential, in the vicinity 
of the gas-solid interface, significantly modifies the elec- 
tron dynamics near the surface, even in the long wave- 
length domain. Here, we have proposed a relatively sim- 
ple approach for the calculation of this spatial variation. 
Future developments will allow for a better estimate of 
its influence on the measurements performed with recent 
experimental techniques like two-photon photoemission 
(TPPE). 
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APPENDIX A: EVALUATION OF THE 
INTEGRAL Ju{z) 

To obtain Ju(z), let us rewrite the vector potential 
A(x,z,t) of Eqs. in section ITTaI as a product of 

two terms related to the variables z, x and t 

po[ A(x,z,t) = pa[ A(z) e <*Z h *-»t) , (Al) 

Then, the interaction Hamiltonian H L (Eq. (jltJafl ) be- 
comes 

h e 



and 



H 



2 mi 



- * ol O(A) eK^x-^t) (A2) 



where 

p O(A) = P A X 

+2M ; 



d d 
2 cos ip— + 2 sin w— + i kf 1 
au ov 

d d p A z 



dz dz 



(A3a) 



s O(A) = 2 *Ay 



■ smp- — h cos<y9— 
ou ov 



(A3b) 



respectively for the p- and s- polarizations. In the equa- 
tions above, the operator poi O(A) is written in the or- 
thogonal crystallographic coordinate system {u,v,z} of 
the surface (see Fig. rotated by an angle p relative 
to the POI coordinate system. The angle p between the 
two coordinate systems can be varied in an experimental 
set-up. For the discrete transitions presented in this pa- 
per, the observables, transition probability and current 
density, arc in principle independent of the angle p. In 
photoemission, one can unravel the local surface symme- 
try (see for example Smith et a/4&). 

The integral Ju(z\ appearing in the laser-matter in- 
teraction matrix elements (|26|) . can now be written in the 
particular case of p- and s- polarizations as 

v JiAz) = ( p Sh(^h+ P ^h^) P shWj (A4a) 



*Ju(z) 



(A4b) 



Taking into account the simplifications related to the 
LWL approximation, one can evaluate these integrals 
along the following lines: The factor Tj" 7 gives a selec- 
tion rule for the angular momenta parallel to the surface 



T}\j = £exp[i(£f 



P 



3=0 



K h )- 



(A5) 



The kP h contribution appears as a consequence of the 
presence of exponential exp(i x) in the spatially vary- 
ing vector potential Eq. (|A1|) . In the LWL approxima- 
tion this exponential is a slowly varying quantity relative 
to the surface reciprocal vectors k u and k v and one can 
use the mean value theorem to factorize it from the in- 
tegrals voi D\j or 9oX S\j. However, because kPj 1 < 
this condition reduces to the standard condition of the 
conservation of the linear momentum fcij ~ fcjj^ parallel to 
the surface. For p- polarization, the explicit expressions 
of the integrals Eq. l|A4a|) are 



*Djj 



ti* i {z')A x (*')Vj{z') dz' (A6) 

(A7) 



2 (cos ip T>™ j + sin ipV}j) 
+i k p x h Sjj 



dA z (z') 
dz 1 



2A(z')^7 
oz' 

Vj(z') dz' 



p 5 



2 Su (cos p + sin p) 



(A8) 
(A9) 



where "Dfj, T>jj and Sjj are function of the reciprocal 
lattice vectors k u and k v . In the equations (|A6|) and 
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one identifies the usual interaction terms present in the 
velocity gauge where the vector potential modulates the 
associated integrals. There are also new terms, called 
surface contributions, originating from the derivative of 
the vector potential, like the last term in the Eq. I|A8|I . 
For s- polarization, the integrals of Eq. (|A4b|l read 



II 
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*Sij(z) 



2 (- sin ip Vjj + cos tp V v u ) (A10) 



tfiz') A„(z')vj(z') dz' (All) 



No surface is term present for this polarization. If one 
retains only the interaction with the near neighbors, the 
inte grals Z >y 7 , Vjj and S u , appearing in Eas.(57)l. l|A"9|l 
and IjAlOfl . present a similar structure. For example, this 
structure can be unravelled by writing explicitly the in- 
tegral Vjj 



u-j-,00,00 



(A12) 



-2cos(fc u a u ){ u Vfj 



+2cos(fc„ al)( u Vfj m + u vff°) 



where the integral "XA^^iA reads 



u 



X }(u-e u a^v~e v a v Q ) (A13) 



Xj{u- t' u a%,v - i' v al) dudv. 

Here the integration domain has been extended over 
three elementary cells in each direction to calculate the 
interaction with the near neighbors. The interaction op- 
erator p0[ O(A) contains symmetric (derivative of the vec- 
tor potential) and antisymmetric contributions relative 
to the inversion of the system of coordinates. This last 
contribution is predominant. Consequently, in Eq. 12211 
the set of basis functions parallel to the surface should 
contain symmetric and antisymmetric basis functions. In 
the simplest case, the functions parallel to the surface 
should contain "s" and "p" like basis functions. 

Finally note that the above simple derivation has 
been inspired by the one given in the appendix F of 
Desjonqueres and SpanjaarcUS where the vector poten- 
tial is independent of the spatial coordinates. 
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